turtles <- read.table('ecol 563/turtles.txt', header=T)
turtles
# reorder the treatment variable
turtles$treat <- factor(turtles$treat, levels=c('fed','fast10','fast20'))
levels(turtles$treat)
# separate groups panel graph
library(lattice)
xyplot(plasma~treat|sex, groups=subject, data=turtles, type='o')
 
library(nlme)
out2 <- lme(plasma~treat*sex, random=~1|subject, data=turtles)
anova(out2)
out1 <- lme(plasma~treat+sex, random=~1|subject, data=turtles)
anova(out1)
summary(out1)
# we saw last time using MCMC that sex should be retained
 
### Graph of marginal model with data ####
 
# panel function version of separate groups plot
xyplot(plasma~treat|sex, groups=subject, data=turtles, type='o', panel=function(x, y, groups,...) {
panel.superpose(x,y,groups,...)})
 
# obtain mean predictions separately by sex
fems <- predict(out1, newdata=data.frame(treat=levels(turtles$treat), sex=rep('female',3)), level=0)
fems
males <- predict(out1, newdata=data.frame(treat=levels(turtles$treat), sex=rep('male',3)), level=0)
males
 
data.frame(treat=levels(turtles$treat), sex=rep('male',3))
 
# add marginal means separately for each sex
xyplot(plasma~treat|sex, groups=subject, data=turtles, type='o', 
panel=function(x, y, subscripts, groups,...) {
 panel.superpose(x, y, subscripts, groups,...)
 sex <- turtles$sex[subscripts][1]
 if(sex=='male') panel.lines(1:3, males, type='o', pch=8, lwd=2, col=1) else
panel.lines(1:3, fems, type='o', pch=8, lwd=2, col=1)
})
 
# change individual trajectory colors to gray70
xyplot(plasma~treat|sex, groups=subject, data=turtles, type='o', 
panel=function(x, y, subscripts, groups,...) {
 panel.superpose(x, y, subscripts, groups, col='grey70',...)
 sex <- turtles$sex[subscripts][1]
 if(sex=='male') panel.lines(1:3, males, type='o', pch=8, lwd=2, col=1) else
panel.lines(1:3, fems, type='o', pch=8, lwd=2, col=1)
})
 
# plot marginal and conditional means together in same plot
# the conditional mean includes the turtle random effects
turtles$pred <- predict(out1)
turtles[1:8,]
ranef(out1)
 
# create a new subject variable that repeats the subject values for each sex
turtles$subject2 <- ifelse(turtles$subject>4, turtles$subject-4, turtles$subject)
 
xyplot(plasma~treat|sex+factor(subject2), data=turtles, type='o')
 
# change order of conditioning levels and the layout
xyplot(plasma~treat|factor(subject2)+sex, data=turtles, type='o', layout=c(4,2))
 
# subject2 needs to be a factor to get useful labels
# when left numeric we get a thermometer label
xyplot(plasma~treat|subject2+sex, data=turtles, type='o',layout=c(4,2))
 
# panel graph of conditional and marginal means
xyplot(plasma~treat|factor(subject2)+sex, data=turtles, type='o', layout=c(4,2), panel=function(x,y,subscripts,...){
panel.xyplot(x,y,col='grey30')
# plot marginal means
sex <- turtles$sex[subscripts][1]
if(sex=='male') panel.lines(1:3, males, lwd=2, col=2) else
panel.lines(1:3, fems, lwd=2, col=2)
# conditional means (include random effects)
panel.lines(x, turtles$pred[subscripts], col=4, lty=2)
})
 
 
xyplot(plasma~factor(treat, levels=levels(treat), labels=c('well\nfed', '10-day\nfast', '20-day\nfast'))| 
factor(subject2)+sex, data=turtles, xlab='Treatment', ylab='Plasma protein (mg/l)', type='o', layout=c(4,2), 
panel=function(x,y,subscripts,...){
panel.xyplot(x,y,col='grey30')
# plot marginal means
sex <- turtles$sex[subscripts][1]
if(sex=='male') panel.lines(1:3, males, lwd=2, col=2) else
panel.lines(1:3, fems, lwd=2, col=2)
# conditional means (include random effects)
panel.lines(x, turtles$pred[subscripts], col=4, lty=2)
})
 
#add legend at bottom of graph
 
xyplot(plasma~factor(treat, levels=levels(treat), labels=c('well\nfed', '10-day\nfast', '20-day\nfast'))| 
factor(subject2)+sex, data=turtles, xlab='Treatment', ylab='Plasma protein (mg/l)', type='o', layout=c(4,2), 
panel=function(x,y,subscripts,...){
panel.xyplot(x,y,col='grey30')
# plot marginal means
sex <- turtles$sex[subscripts][1]
if(sex=='male') panel.lines(1:3, males, lwd=2, col=2) else
panel.lines(1:3, fems, lwd=2, col=2)
# conditional means (include random effects)
panel.lines(x, turtles$pred[subscripts], col=4, lty=2)
}, key=list(x=.8, y=-.15, corner=c(0,0), lines=list(col=c('grey30',4,2), pch=c(1,NA,NA), lty=c(1,2,1), lwd=c(1,1,2),
size=2, type=c('p','l','l'), divide=2, cex=.85), text=list(c('observed', 'conditional mean', 'marginal mean'), cex=.75), 
border=T))
 
# save graph
xyplot(plasma~factor(treat, levels=levels(treat), labels=c('well\nfed', '10-day\nfast', '20-day\nfast'))| 
factor(subject2)+sex, data=turtles, xlab='Treatment', ylab='Plasma protein (mg/l)', type='o',layout=c(4,2), 
panel=function(x,y,subscripts,...){
panel.xyplot(x,y,col='grey30')
# plot marginal means
sex <- turtles$sex[subscripts][1]
if(sex=='male') panel.lines(1:3, males, lwd=2, col=2) else
panel.lines(1:3, fems, lwd=2, col=2)
# conditional means (include random effects)
panel.lines(x, turtles$pred[subscripts], col=4, lty=2)
}, key=list(x=.8, y=-.15, corner=c(0,0), lines=list(col=c('grey30',4,2), pch=c(1,NA,NA), lty=c(1,2,1), lwd=c(1,1,2),
size=2, type=c('p','l','l'), divide=2, cex=.85), text=list(c('observed','conditional mean', 'marginal mean'), cex=.75), 
border=T))->mygraph
 
# print graph with side and top panels
library(latticeExtra)
useOuterStrips(mygraph)